!!! This module includes common functions for setting spline
MODULE comf_setspline
  use nr, only : spline  !using Numerical Recipe (NR) subroutines
  IMPLICIT NONE       
CONTAINS
    
  !==========================================================
  ! Initialize the cubic spline, one dimension over xgrid_
  ! linear interpolate the slope at both ends.
  !==========================================================
  subroutine set_spline(xgrid_,y_,spgrid_)
     IMPLICIT NONE
     real(8), dimension(:), INTENT(IN) :: y_, xgrid_
     real(8), dimension(:), INTENT(OUT) :: spgrid_
     real(8) :: xp(3), yp(3), yp1,ypn,rtemp
     integer :: iSIZE
     xp(1) = y_(1)
     xp(2) = y_(2)
     xp(3) = y_(3)
     yp(2) = (xp(2)-xp(1))/(xgrid_(2)-xgrid_(1))
     yp(3) = (xp(3)-xp(2))/(xgrid_(3)-xgrid_(2))
     rtemp = (xgrid_(1)-xgrid_(2)) / (xgrid_(3)-xgrid_(2))
     yp1 = rtemp * yp(3) + (1.0d0-rtemp) * yp(2)
     
     iSIZE = size(y_)
     xp(1) = y_(iSIZE-2)  !vV%a(git,giA,N_S-2,giAIME,giHEALTH,giRECSS)
     xp(2) = y_(iSIZE-1)
     xp(3) = y_(iSIZE)
     
     yp(1) = (xp(2)-xp(1))/(xgrid_(iSIZE-1)-xgrid_(iSIZE-2))
     yp(2) = (xp(3)-xp(2))/(xgrid_(iSIZE)-xgrid_(iSIZE-1))
     rtemp = (xgrid_(iSIZE)-xgrid_(iSIZE-2)) / (xgrid_(iSIZE-1)-xgrid_(iSIZE-2))
     ypn = rtemp * yp(2) + (1.0d0-rtemp) * yp(1)
     call spline(xgrid_,y_,yp1,ypn,spgrid_)
 end subroutine
 
   !==========================================================
  ! Initialize the cubic spline, one dimension over xgrid_
  ! shrink one grid at each end to get slopes.
  !==========================================================
  subroutine set_spline_shrinked(xgrid_,y_,spgrid_)
     IMPLICIT NONE
     real(8), dimension(:), INTENT(IN) :: y_, xgrid_
     real(8), dimension(:), INTENT(OUT) :: spgrid_
     real(8) :: yp1,ypn
     integer :: iSIZE
     
     iSIZE = size(y_)
     yp1 = (y_(2) - y_(1)) / (xgrid_(2) - xgrid_(1))
     ypn = (y_(iSIZE) - y_(iSIZE-1)) / (xgrid_(iSIZE) - xgrid_(iSIZE-1))
     
     call spline(xgrid_(2:(iSIZE-1)),y_(2:(iSIZE-1)),yp1,ypn,spgrid_)
 end subroutine
   
END MODULE comf_setspline
